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^ Abstract 

l/~j . The solution to the BFKL equation grows like a power of center of mass energy, 

s, violating unitarity conditions at high energies. The growth of the cross section 
can be tamed by taking into account multiple pomeron exchanges. This is known as 

■^j- ■ saturation and it is expressed in the Balitsky-Kovchegov equation, [Sj. Conservation 

of energy should also slow down the growth of the cross section, and our aim in this 

,-SZh ' work is to study the effects of enforcing energy conservation in DIS-events. Using 

i 1, the dipole picture, onium-onium collisions can be viewed, in the large N c and the 

Qh. leading logarithmic limits, as the scattering of a collection of color dipoles. We 

construct a Monte Carlo program, based on Mueller's model, El IU IE] > and using 

energy conservation, to study onium-onium and onium-nucleus collisions. Dipole 



fusion processes will be important at high gluon densities and should also slow down 
the growth of the cross section. We propose an expression for the fusion factor and 



use it in our Monte Carlo to study its effects. 
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Figure 1: Deep Inelastic scattering between an electron and a proton. 

1 Introduction 

It remains that, from the same principles, I now demonstrate the frame of 
the System of the World. 
Isaac Newton 
Book III, Philosophiae Naturalis Principia Mathematica, 1687 

Deep Inelastic Scattering (DIS) is a scattering process where the target particle disin- 
tegrates and what comes out is very different from what came in. If on the other hand 
the system after the collision looks exactly like the system before the collision we call the 
process elastic. A very common process to study is the collision between an electron and 
a proton. This process is usually deeply inelastic at high energies, which is due to the 
fact that the proton is not an elementary particle but is instead made up of even smaller 
particles, the quarks. This was observed in 1968, and it was found that the proton consist 
of three quarks, two up-quarks and one down-quark. However, the picture is not as simple 
as that. Since the quarks interact via the strong force which is mediated by the gluons, 
there are also a lot of gluons in the proton. These gluons can create quark-antiquark pairs 
or split into new gluons, and thus the proton structure is much more complicated. It was 
indeed observed in these collisions that almost half of the proton's momentum is carried 
by particles that don't interact electromagnetically, namely the gluons. Gluons and quarks 
are collectively called partons. 

Let us consider a typical electron-proton scattering event, see figure 1. We have an 
incoming electron with four-momentum p and an incoming proton with four-momentum 
P interacting via the exchange of a virtual photon with virtuality Q 2 = —q 2 . 

With P' we denote the four-momentum of the parton which absorbs the virtual photon. 
The simplest possible process is that it stays on the mass-shell after the absorption, with 
a momentum P' + q. We define the variable x, as the fraction of the proton's momentum 
carried by the parton, which is struck by the virtual photon; thus x = P'/P. Since the 
mass of the proton, and naturally its partons, can be neglected at the energy scales we are 



concerned about, we get (P' + q) 2 = and P' 2 + IP' ■ q + q 2 = IxP ■ q — Q 2 = 0, and hence 
x = j^-. This variable x is called the Bjorken scaling variable, usually denoted by xbj- 



2 Theory 

We define, the parton distribution function f(x), as the probability for finding a parton 
carrying a momentum fraction x of the proton. Usually we use subscripts to denote which 
parton we are considering, for example with f q (x) we denote the probability for finding a 
quark with momentum fraction x. To specify the flavor of the quark we use the subscripts 
u,d etc. The parton model predicts that these distribution functions depend only on the 
variable x. This behavior is known as scaling, which is also why x is called the Bjorken 
scaling variable. If a parton has a certain x it could either have this momentum fraction 
from the beginning, or it could initially have a larger momentum fraction and reduce it 
to x through emission of gluons . The available phase space for radiating gluons is larger 
when Q 2 is large. Therefore the distribution functions will show a greater probability to 
have smaller x, and smaller probability for carrying large x for larger Q 2 . This implies 
that the distribution functions are not only functions of x, but they are functions of x and 
Q 2 , f(x) = f(x,Q 2 ). This is a central feature of QCD, and it is called scaling violation. 
It is possible to derive evolution equations for the distribution functions in the variable 
Q (logQ is often used instead) if one considers the various processes in the proton, such 
as quarks radiating gluons or gluons splitting into gluons or quark- ant iquark pairs. Thus, 
if one knows the distribution function for a certain value of Q, the value for any other 
Q can be calculated using these equations, which in the case of QCD are known as the 
Altarelli-Parisi Equations. These equations can be found in any book on Quantum Field 
Theory. 

The distribution functions are important because the cross sections for different events 
(both in QED and QCD) can be expressed in these function. Thus, for the process we 
considered above, an electron hitting a proton, the cross section can be written as 

° = Yl / dx \.fi( x ' Q 2 ) a ( e ~Q -»• e ~<?) + fq(x, Q 2 )o-(e'q -»• e~q)\ (1) 

M 

Here the cross sections in the integrand are the cross sections for the sub-collisions between 
the electron and the partons, and the sum runs over the quark and the antiquark flavors. 
From now on we will use capital letters to denote the densities in the variable lnl/x, 
F(x, Q 2 ) for fermions and G(x, Q 2 ) for gluons. 

2.1 DGLAP, BFKL and the LDC Models 

2.1.1 DGLAP evolution 

Consider now the emission of many gluons before the parton is hit by the photon. This is 
represented by a fan diagram, see figure 2. We define Zi as, Zi = k+i/k+i-i, hence 1 — Zi 








Figure 2: A fan diagram for a DIS event. The virtual propagators are denoted by ki and 
the quasi-real gluons from the initial-state radiation are denoted by q$ . 



is the fraction of positive light cone momentum, p + — Po+Pi, carried away by the quasi- 
real gluon, qi. This implies that the Bjorken variable, x, is given by x = YTt=i z i x o- I n 
the double leading log approximation the gluon ladder is strongly ordered in transverse 
momentum 

qli<q 2 ±2<---qin<Q 2 (2) 

and in s, 2j C 1. The probability, in the leading log approximation, to emit a gluon is 
given by 



Aa s dzi dq\- 



(3) 



3tt 1-Zi q\i 

This factor must however be modified because the ordering in (J2J) implies that no gluon 
is emitted with qj_ between q± t i-i and q± t i. Thus we must also take into account the 
probability that no emission has occurred between these two values. This correction factor 
is called a Sudakov form factor, and it is, in this case, given by 



S{Q±,ii <?_L,j-l, 
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(4) 



where the e is introduced as a cutoff since the integral diverges for z = 1. Generally, if 
we have an ordering in a variable t and the probability, P(t), for something to occur at a 
certain t value, the Sudakov form factor is given by 



S(t u h 



exp 



t-2 



dtP(t) 



(5) 



Here (ti,^) is the interval where an emission was allowed but did not occur. Summing 
over the number, n, of links in the chain, the following expression for the structure function 
can be derived (for a quark chain) 

F ~ £ fl f / ^T ^T^fc' ^-i)%^-^-i)l S(x-f[ z jX o)0(<?-<l in ) (6) 



n i=l 



Subleading corrections in lnl/rr are taken into account by replacing the 1/(1 — z t ) term 
by the full splitting function, P(z). Including spin effects gives 

i U + ." (7) 



1-z 2 1- z 

By taking the derivative of eq. © with respect to InQ 2 we obtain the DGLAP (Dokshitzer- 
Gribov-Lipatov-Altarelli-Parisi) evolution equation 



dF(x,Q 2 ) 4a 



<91nQ 2 3tt 



j^H'W.cm*-*!) 



Here we have replaced P with P (we describe P below), and using P one can let e — > 
0. This can be done since the divergence can be treated by defining a distribution that 
can be integrated by subtracting a delta function from the singular term. We define the 
distribution, , ± K , such that for any given smooth function f(z) we have 

X d2 (r^y: - X d2 d- Z ) (9) 

and then we replace the splitting function P with P given by 

MfT^H* 1 -*) (10) 

In fact, it is P and not P that appears in the AP equations. The delta function comes 
from the fact that we have to take into account that a parent quark splitting into a quark 
and a gluon disappears. Thus, we get contributions from the new particles but the one in 
the initial state no longer exists, and this must be taken into account. The coefficient in 
front of the delta function, 3/2, is found by the requirement that the number of quarks is 
conserved, which means that the integral of P over z must be zero. In the case of gluons 
(i.e a process where a gluon splits into two gluons) the number of gluons is not conserved 
but the total energy is, therefore, in this case, the integral of zP over z is zero. 

2.1.2 BFKL region 

Assume that Q 2 is only moderately large and that x is kept small. The splitting function 
for the process where a gluon splits into two gluons, P gg , contains a pole, 1/z. This implies 



that for small x and moderately large Q 2 , this process is very probable and will dominate 
over splitting events involving quarks. Therefore we now consider gluonic chains where a 
parton emits a gluon which then develops a chain of gluons by successive splittings into new 
gluons. In this region, non-ordered ladders, though suppressed, will be more important 
and they have to be taken into account. This means that the transverse momenta, k±, can 
now also decrease. Considering this new possibility, Balistky, Fadin, Kuraev and Lipatov 
derived a more complex evolution equation, which is known as the BFKL equation. The 
evolution in x is given by 

dF(x,kl) f ,,, 2 



<91n(l/x) 



I dk'lK(k 2 ± ,k'l)F(x,k'l) (11; 



Here we have introduced the unitegrated distribution function, T. There are various 
definitions of J 7 , and in the LDC model (which is discussed in the next section) it satisfies 
the following relation 



—j- {x , k±)+ j—j- {x — k±) - 



F(x,Q 2 )= I ^tH^kl)+ ^tH^kl)^ (12) 



The kernel of the integral, K, in equation (jllj) . is known as the BFKL kernel. The 
important thing here is to find the eigenvalues of this kernel. To leading order in ln(l/x), 
and for a fixed coupling constant as, the largest eigenvalue is given by A = 12aslog(2)/7r. 
It can then be shown that the gluon distribution function is given by 

G ^(Aiog(i/,)r =i _ A 

z — ' n\ 

n 

We will come back to the BFKL equation, and the so called BFKL pomeron, later when 
we discuss saturation. For now, one can note that BFKL and DGLAP are applicable in 
different regions. In the next section we will describe two models which interpolate between 
these two different regions, thus containing both the DGLAP and the BFKL physics. 

2.1.3 CCFM and the LDC Models 

The first such model was proposed by Catani, Ciafaloni, Fiorani and Marchesini and is thus 
called the CCFM model. We will discuss purely gluonic chains and the last gluon in the 
chain is assumed to interact with a color neutral probe. One can usually order emissions 
in different variables, such as color, rapidity, transverse momentum and so on. Due to soft 
color coherence |T3j the ordering in rapidity and in color are equivalent however. This is 
used in the CCFM model where the initial state radiation is made up of gluons which are 
not followed, in rapidity, by a gluon with more energy, or equivalently positive light-cone 
momentum, q + = q + q\. 



With the selection of the initial state radiation described above, it was shown by the 
authors that the unintegrated distribution function is given by 

Q(x,kl,q) ~ J2i[\[®^(- A ne(z l ,kl l ,q l ) + -L-)A s x (14) 

„ i=1 U <2±,i \ Z i L Z iJ 

n 

Q(q~i ~ Qi-iZi-i) S(x - J J Zj)9(q - q n z n )5(k]_ - k 2 ±n ) 

J 3 

where q\ is defined by g, = q± j i/(l — z.j). A ne stands for a non-eikonal form factor and A5 
is the Sudakov form factor. They are given by 

!, M a ( - f dq h dz 

- ln — ; A 5 = exp I -a / -J— 



A ne = exp( -aln-ln— ^ ); A 5 = exp( -a / -3-^- -0 or der ) (15) 



Here Qorder specifies the kinematic region allowed by the ordering constraints, and a = 
3a s /ir. In the CCFM model, only the singular terms, 1/z and 1/(1 — z), in the splitting 
functions are included. The 1/z pole is important in the BFKL region and the 1/(1 — z) 
pole and the Sudakov factor are more important in the DGLAP region. Since the non- 
singular terms in the splitting function are not included, DGLAP is not fully reproduced 
in the large Q 2 and the large x limits. 

The LDC (Linked Dipole Chain) model, [TJ, is a reformulation and a generalization 
of the CCFM model. The main idea in this model is that certain gluons in the initial 
state radiation of the CCFM model, namely those gluons which satisfy the constraint 
q±i < min(g_ L (j_ 1 ), q±(i+i)), can be treated as final state emissions from dipoles, which are 
formed by the remaining initial state gluons. These remaining gluons, the initial state 
gluons in the LDC model, will then satisfy q±i ~ max(/cj_j, fcj_(j_i)). They will be ordered 
in both positive and negative light-cone momenta, q- — qo — qi- A single chain in the LDC 
model corresponds to several chains in CCFM, and it was shown in pQ that the non-eikonal 
factors in CCFM add up to one, when one considers all contributions from these chains. 
The ordering, in both p + and p_, and the absence of the non-eikonal form factors, implies 
that LDC is forward-backward symmetric. If one starts from the probe end, which has 
large g_, and follows the chain towards the proton end, which has large P + , the situation 
will look the same, the minus and plus components will just change roles and we will have 
exactly the same ordering and the same distribution function. In the LDC model the 
distribution function satisfies the following evolution equation (see eq.43 in [T]) 



<9jF(x, k 2 , ) 

— fa a 

<91og(l/x) 



r k ± Hk' 2 , r rih' 2 h 2 h' 2 , 

J k' 2 { ' ±J y fe2 k' 2 k' 2 { kV L 



2 U2 u'2 

(10) 



Here k\ is the transverse momentum of the last gluon in the chain while k 2 is the momen- 
tum in the step before. The first term here is the DGLAP contribution where the transverse 
momentum increases in the last step, while the second term contains the probability that 
the transverse momentum decreases in the last step. For a more detailed discussion on 
these models, see [TUj . 





Figure 3: A quark and an antiquark emitting a gluon. In the first picture the red quark 
radiates a red-antiblue gluon and changes its color to blue and in the second diagram the 
antiblue antiquark emits a antiblue-red gluon and changes its anticolor to antired. 

2.2 Mueller's Dipole Formulation 

The idea of dipoles as emitters of radiation is quite natural in QCD. Consider the two 
diagrams in figure 3. The two diagrams will interfere and the two final states are identical. 
Thus, it is not possible to tell whether the quark or the antiquark emitted the gluon. 
Therefore it makes sense to regard the quark- antiquark pair as a color dipole, radiating 
coherently. In the region between the quark and the gluon, the blue and antiblue color 
separation will give rise to the emission of softer gluons. But in directions further away the 
emissions from blue and antiblue will interfere destructively. Therefore, in these directions, 
the emission of softer gluons will correspond to a red antired color dipole. In the restframe 
of the quark and the gluon however, emission of softer gluons corresponds to a blue antiblue 
color dipole. Therefore the emission of softer gluons corresponds to two independent color 
dipoles. This generalizes and a system of a quark, an antiquark and n — 1 gluons will 
correspond, in the large N c limit, to a system of n independent color dipoles. Each time 
a gluon is radiated a dipole splits into two new dipoles and the process continues. This is 
used in the Dipole Cascade Model, J21EIJ; which is formulated in momentum space. 

A large part of this thesis is based on Mueller's Dipole model J2J El IH E] and we will 
be using it in our MC program. In Mueller's model we start with a heavy quark-antiquark 
pair, an onium state, and consider successive emissions of softer and softer gluons. The 
main idea here is also based on the observation that, in the large N c limit and in the leading 
logarithmic approximation the complete, multi-gluon, wavefunction of an onium state can 
be regarded as a collection of dipoles. This is so since, as we saw above, each gluon will 
act like a quark-antiquark pair. Mueller's approach is formulated in transverse position, 
x, instead of transverse momentum, k. These two variables are naturally connected by a 
Fourier transform 

,(o)f_ v ^_ f d2k i B ai-xi./,C)i 
(2nf 

Here ip is the light-cone wavefunction (the zero on ip indicating that there are gluons 
present) and a and (3 are heavy quark and antiquark spinor indices, xo and xi are the 
transverse coordinates of the quark and the antiquark respectively. Finally, Y is the total 
rapidity range determined by the center of mass energy, Y ~lns/M 2 , M being the mass of 



W*i> y ) = / 7^ eikl ' Xl CJ(ki,n (17) 




Figure 4: Creation of new dipoles by means of gluon emission. The original quark is 
denoted by while the antiquark is denoted by 1. The original dipole "decays" into two 
new dipoles when the gluon labeled by 2 is emitted. The process continues with the emissions 
of gluons 3 and 4 etc. 

the onium. We also define Xio = xi — xo, which means that the transverse size of the initial 
dipole is, x\q = |x 10 |. The reason we choose to work with transverse coordinates instead 
of transverse momenta is that one can assume the partons to be frozen during emissions 
and the evolution process, which makes the analysis much easier. This is justified since 
the time scales for the processes that we consider are very short. Denoting the square of 
the wave function with </>, Mueller showed that 
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where x 2 is the transverse position of the emitted gluon and (p^ and (p^ are the squared 
wave functions for the onia state without and with a gluon respectively. The factor 

o 

/2 x 

« X 2ZT 



- plays the same role as the splitting functions we hade earlier. Observe that 
the X2 integral diverges for small £20 or x i2- This means that the probability to emit a 
dipole goes to infinity as the dipole size approaches zero. To deal with this a cutoff, p, is 
introduced, which later can be sent to zero and will not appear in any physical quantities. 
Generally, if we have n dipoles and we consider the emission of the nth gluon, we will get 
a contribution 



E/d 2 * 



.?;- 



n 1 1 



(19) 



where ij runs over all the indices for which there exist a dipole. 



The complete formula for the wave function with arbitrarily many gluons is derived 
from a generating functional. Mueller derived a non-linear integral equation for this gen- 
erating functional. Taking into account not only real gluon emissions, but also the leading 
logarithmic virtual corrections, the formula reads 



Z(b ,Xio,Y,u) = u(h ,xL 10 )exp 

<-Y 



-*Ht) y . 



+ 



+ 7T I d V ex P 
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f x 2 1 1 

x / rf 2 x 2 2 10 2 Z(b + -x 12 , x 20 , y, w)Z(b - -x 20 , xi 2 , y, u) 

J p X 2§ X Y2 *• *• 

Here the p under the second integral indicates that we are integrating over a region satis- 
fying, rr 20 ^ p and x 2 \ ^ p. The vector b is the position of the center of the initial dipole 
in transverse coordinates. The two exponentials have the meaning of Sudakov factors, 
where in the first one we consider the probability that nothing happens within the allowed 
rapidity interval, which means that there are no gluons present. The second exponential 
contains the probability that nothing happens after the last gluon is emitted at rapidity y. 

2.3 Saturation and The Balitsky-Kovchegov Equation 

Consider the scattering of two onium states against each other. Using the dipole formula- 
tion described above, this process can be viewed as the scattering of a collection of color 
dipoles against each other. The single pomeron exchange then corresponds to the approx- 
imation in which two dipoles (from the two different states) interact via the exchange of a 
color neutral gluon pair, independently of the rest of the dipoles. This pomeron is known 
as the BFKL pomeron. The total cross section would then be given by the number of 
dipoles in one onium state times the number of dipoles in the other onium state times the 
cross section for dipole-dipole scattering due to the two gluon exchange. However, as the 
center of mass energy increases the number of dipoles will increase, and when the number 
becomes sufficiently large the total cross section will violate the unitarity bounds. The 
reason is that at higher energies, second order processes must be taken into account, and 
these will slow down the growth of the total cross section. The higher order scatterings 
are called multiple pomeron exchanges. 

Consider the scattering of a photon off a nucleus. We assume that the photon splits 
into a quark- ant iquark pair, on which we apply the dipole formulation, that is we work in 
the large N c and the leading logarithmic limits. This means that the photon will develop 
a cascade of color dipoles, which interact with the nucleus. The nucleus is assumed to be 
at rest. We assume that each dipole interacts with the nucleus by a two gluon exchange, 
independent from the rest of the dipoles. The nucleus is approximated by a cylinder which 
makes the calculations easier. The multiple pomeron exchanges will correspond to the 
simultaneous interactions of multiple dipoles with the nucleus, which means that the fc:th 
pomeron exchange corresponds to k dipoles simultaneously interacting with the nucleus by 

10 



exchanging k color neutral gluon pairs. With these considerations and summing pomeron 
exchanges up to all orders, Kovchegov, in j^j, derived the following non-linear differential- 
integral equation for the forward scattering amplitude off the nucleus: 

dN(-x 10 ,b ,Y) a [ 2 Ao f AT , v , 1 v x, 

— w — = 2^y/ X2 ^f;( v iV(X20 ' bo + 2 Xl2 ' r) + 

+iV(x 12 , b - ix 20 , Y) - iV(x 10 , b , Y) - (21) 

1 1 

-N(x 20 , b + -x 12 , y)JV(xi2, b - -x 20 , Y) 

Here the imaginary part of N gives the total cross section. This equation is known as 
the Balistky-Kovchegov (BK) equation (a very similar equation was proposed by Balitsky 
in |llj). and as we mentioned above, it contains all the multiple pomeron exchanges. The 
linear part of the equation is the BFKL part which dominates for small values of N and 
which causes the cross section to rise exponentially. The non-linear part comes from the 
multiple pomeron exchanges and causes the amplitude to saturate for values near 1. The 
equation has a nice physical interpretation, the first two terms on the right hand side are 
the contributions to the scattering amplitude from the creation of the new dipoles, 12 and 
20, while the third term comes from the fact that the original dipole, 10, no longer exists. 
The reader might guess that the non-linear term is the contribution from the recombination 
of the two dipoles, 12 and 20, into the original dipole 10, but as we said above the non-linear 
term comes from the multiple pomeron exchanges. The Balitsky-Kovchegov equation does 
not contain the process of two dipoles recombining into one. This fusion process should 
also contribute to saturation. For large gluon densities it should be non-negligible, but it 
is not an easy task to estimate this effect. In fact, in the derivation of (|21|). Kovchegov 
used (|20|) . which does not contain the fusion process of dipoles. We will come back to 
dipole fusion processes on section 5. 

3 MC Simulation of Onium Evolution With Energy 
Conservation 

3.1 The Main Idea 

As we mentioned earlier, the probability to emit small dipoles becomes infinitely large as 
the dipole size approaches zero. To deal with this we introduced a cutoff, p. This cutoff is 
an ultraviolet cutoff since a small dipole corresponds to two hard gluons. This follows from 
the relation - ~ px, where r is the transverse size of the dipole and p± is the transerverse 
momentum. We will now develop an onium state using an ordering in rapidity, y, such that 
the rapidity increases for each emission. We find it convenient to define y =ln— , hence 
p + = p^e~ y . We see that, for fixed y, p + — > 00 as r — > 0. This means that, if we allow 
very small r we will violate energy conservation since the energy of the gluon will become 

11 




Figure 5: The photon develops a cascade of dipoles which interact with nucleus by the 
exchange of color neutral gluon pairs. Each double line represents a gluon which behaves 
like a quark- antiquark pair in the large N c limit. 
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larger than the energy of the mother dipole. Hence, we cannot allow too small dipole sizes. 
The positive light-cone momentum is essentially equal to the energy, and from now on we 
will refer to it as energy. Let p + be the energy of the parent dipole and let the gluon be 
emitted at a rapidity y. Denote by p +g the energy of the gluon. We then have 

1 e~ y 

P+q ^ P+ =^ P±o e y = -e y ^ p + , hence r ^ (22) 

r p + 

This means that we get a very natural cutoff for small dipoles, given by, p(y) = e~ y /p + . 
To make the program right-left symmetric, which means that the theory looks the same 
independent of which way we evolve, we must also demand p- conservation, which is 
described later. These observations are the main motivations for this thesis. We will now 
describe the Monte Carlo program, including energy conservation, developed to study the 
deeply inelastic scattering between a photon and a nucleus, and between two onium states. 

3.2 Recoil Formulas 

Since we demand energy to be conserved, we must also take into account the full recoil 
effects. A single gluon will, in the dipole model, be a part of two independent dipoles and 
the question is how the two dipoles share the recoil from the emission. We solve this by 
assuming that when a dipole emits a gluon, the initial gluons, which make up the parent 
dipole, may contribute with all their energy. This means that the next time a neighboring 
dipole emits a gluon, the available energy is reduced because one of its gluons has lost 
energy from an earlier emission. How much energy does each parton in the dipole give to 
the emitted gluon? Obviously the answer should be such that if the new gluon comes very 
close to one of the partons, that parton should take most of the recoil, while the other one 
should be less affected. There are, however, different ways to share the amount of energy 
between the partons. 

Consider the emission of gluon n from the dipole ij. Let r« n and r n j denote the distance 
between the new gluon and parton i and parton j respectively. Denote the recoil by p" , 
while p' + is the amount of energy left after the emission, which implies that p + = p" +p' + . 
We then make the ansatz that 

P+i = ; P+n and p +j = ■ p +n (23) 

1 jn ~r ' m ' jn t ' in 

Thus we assume that p\Jp" + j = rj n /r in . We have also studied alternative ways to share 
the recoils but the result does not depend sensitively on the exact formula. When we have 
generated an emission we always check that p" ^ p + , which obviously must be satisfied. 

To make the analysis simple, transverse momentum, pj_, is only approximately con- 
served in the program. Exact conservation of pj_ would be possible but is not essential 
for the purpose of this study. When a new gluon is emitted its transverse momentum is 
decided by 

p Ln = max I — , — ) (24) 

\'in 'jn/ 

13 



The recoils on the emitters are given by 

If p ±l < — then p' ±i = — (25) 

Tin '"in 

lfp±i ^ — thenp',- = — 

i ]n i jn 

Otherwise we simply set p± = p' ± , i.e we don't change the transverse momenta. This simply 
means that the transverse momentum is decided only by the largest contribution. Hence 
the transverse momentum of a parton will always be determined by the shortest distance 
to another parton, with which it has formed a dipole. Since we change p + and p± we must, 
to be consistent, also change the rapidity y. This is obviously done by setting 

y' = \J±>y (26) 

P+ 

The inequality is due to the fact that p' ± ^ p± and p' + ^ p + . Therefore the rapidity 
will always increase due to a recoil. There is however a problem with this; as we order 
the emissions in rapidity the last emitted gluon will always have the largest rapidity, and 
we would like it to remain so. Since the rapidity of the emitters increase due to recoil, 
there is the possibility that one of the emitters, or both of them, after the recoil ends 
up with a rapidity larger than the rapidity of the emitted gluon. The next time we emit 
a gluon, the rapidity of that emitter would then determine the region in which the new 
emission can occur, in a sense violating the ordering we had. To prevent this we require 
that y'^y'j ^ y n ] if this constraint is not satisfied a new gluon is generated. Without this 
constraint there would also be a lot of gluons in the final state which have rapidities larger 
than the maximum rapidity, due to the recoils. 

3.3 Upper Bound for the Dipole Size 

The p + conservation, as we mentioned, will give a cutoff for small dipoles, limiting the 
number of small dipoles. We can, however, also get extremely large dipoles, especially when 
the available energy is low and the rapidity is not very large. This follows immediately 
from the fact that p(y) = e~ y /p + gets very large for limited y and small p + . We can also 
take extremely small steps in rapidity because the phase space available permits this, see 
figure El When we are about to emit a new gluon we let the computer generate an emission 
for each existing dipole, and then we choose the emission which has the lowest y. Because 
of what we said above, the step sizes in y can, and indeed will, be extremely small which 
result in a huge number of dipoles and very long running times. This can be fixed by 
imposing an ordering in negative light cone-momentum, p_, such that 

p- n ^ max(p^i,p-j) (27) 

There is also a physical argument for doing this. As we described in section 3.1 we 
want the cascade to be right-left symmetric, and by imposing p_ ordering we can achieve 
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Figure 6: The triangular phase space for emissions. The two circles indicate the emitting 
partons, the vertical line through circle j indicates the rapidity of j while the on second 
line, p- is equal to p^j Emissions between the two lines through j are forbidden due to 
P- ordering, which limits the allowed rapidity interval. The line through i indicates its p + 
value. 

this. This ordering will make very small steps in y less probable, and it will also limit the 
radius of a dipole from above. To see this let the larger negative light-cone momenta be 

p_ which implies 

1 e y 

^ P- =>- p±e y = —e y ^ p- hence r ^ — (28) 



P- 



-e y ^ p. 

r 



hence r ^ — 
P- 



We are limiting the small dipoles with the energy available in the emitting dipole, so 
the smallest dipole that can be created is the same for both partons of the mother dipole. 
A more realistic way might be to set two different lower limits for the two partons. If one 
of them has greater energy than the other, it should be possible for the emitted gluon to 
form a smaller dipole with this parton than with the other one. Thus, the energy of the 
gluon should be bounded by the energy of the two partons separately, and not by their 
sum. We would then have 



e- y , e~ v 

^ and r jn ^ 

P+i ' P+j 



(29) 



This is however already satisfied in our program due to our constraint y[ ^ y n , which is 
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necessary for conserving the y ordering. The constraint 



-ln- 



y[ ^ y n implies e Vi ^ e y ". The relation e y% = e p ^ 
1 



— — then implies 
P±i 



-yn 



-yn 






> 



P±i P+i P+i 
1 1 



which together with the relation 



Ai 



max\ p±i 



gives r in ^ 



P+i 



^ n 



(30) 



3.4 Formulas for Generating y and r 

We have above described the main ideas of our program. In this section we will write down 
the explicit equations used to generate y and r values. As we mentioned in section 2.2, 
the r x and r y factors in (J18j) have the meaning of a splitting function. Taking into account 
also the Sudakov factor, the total formula for emission is given by 



Z7T Jp{y") r in r jn 



■exp 



Z7T Jy' Jp(y") 



d 2 r r . 



in jn 



(31) 



The Sudakov factors in section 2.2 were trivial, since we had no y dependence on the 
r integral. Demanding energy conservation implies, however, that the lower limit of the 
integral depends on y and the Sudakov factor will therefore be more complex. To make 
the calculations simpler, we set rj = (0,0) and Yj = (1,0), so that r^ = 1. This we can 
always do by applying a translation of — r-j, then rotating the vector r'j = r'^ onto the 
x-axis and finally by scaling with r^. After we have generated values for r we transform 
the system back to the original one. To generate the y and r values we use the well known 
veto algorithm, and the y values are given by the distribution 



y 



- ln \/i + Y + V( ln V 1 + 4 + 2/,)2 



IilR 

4a 



(32) 



where R is a random number, y' is the maximum rapidity before the emission, and p + 
is the total energy of the emitting dipole. The values for r x and r y are chosen from the 
distribution 



P 2 (l + i^ 2 



V 



\ i + V_v ( i + 2 i,)« 2 



P 2 (1 + 4^) R2 



\jl +V _ V(1 + _^)R 2 
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cos 2tcRi 



sin 2ttRi 



(33) 



where p = p(y) and R± and R<i are two, independent, random numbers. Note that we use 
the generated y- value for generating r x and r y . Finally, we accept the generated values 
with the probability 

ln(l + P ^e 2y ){{r x - if + r \ + 0.25)(r" + r 2 y + 0.25) 

(ln(l + P -±) + 2y) f((r x - If + rl)((r x - 1)2 + r* + 0.25) + (7* + r^ + 0.25)(r2 + r %)\ 

(34) 
This formula is derived in the appendix. 

3.5 Summary of the Program 

In this section we will summarize the important points from the earlier sections about the 
MC program. The algorithm is 

1. For each dipole, pick a y value from (|32|) and r x and r y from (|33|) . 

2. For these values calculate the weight in (|34|) . If (|34j) < R go back to 1 and start 
over. Otherwise go to 3. 

3. If p +i ^ p+j or P+i ^ P+j g° back to 1 and start over. Otherwise go to 4. 

4. Check (|27|). If the inequality is not satisfied, then go back to 1 and start over. 
Otherwise go to 5. 

5. Check if y[, y'- ^ y. If the inequality is not satisfied, then go back to 1 and start over. 
Otherwise generate a gluon. 

6. Set the new values for p+ij, Pu,j and ytj. 

7. Select the dipole which gives the lowest y. 

8. If y < y max let the selected dipole emit the generated gluon. Otherwise stop the 
process. 

9. Delete all other generated gluons and continue. 

4 Applications of the MC Program 

In this section we consider the applications of the program. We use the program to calculate 
cross- sect ions for two processes, onium-nucleus scattering and onium-onium scattering. 
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4.1 Onium-Onium Scattering 

We start with onium-onium collisions. There has been other studies of this process, in 
particular there is another MC program called OEDIPUS, developed by Salam, [7J, which 
is also based on Mueller's Dipole formulation, but without energy conservation. To study 
this process we generate two onium states independent of each other, one developed up 
to a rapidity y and the other developed up to Y — y, and after that we let them collide. 
The collisions will occur between the dipoles from the two onium states and to calculate 
the cross section we need an expression for the S matrix. It is here in order to present the 
theoretical background for how we arrive at an expression for the S matrix. The amplitude 
for one pomeron exchange between two colliding onia is, in the BFKL approximation, given 
by 

f d 2 cd 2 c' 
F (1) (r l5 r 2 , b, Y) = - / ^ a d 2 rd 2 r'f(r - r, c, c')n(c, n, r, y)n(c', r 2 , \r' -b\,Y - y) 

(35) 
Here r± and r 2 are the sizes for the two initial dipoles, c is the size of a dipole in the 
right moving onium while d is the size of a dipole in the left moving onium and r and 
r' are their positions relative the initial dipoles respectively. n{c,r\,r,y) is the density of 
dipoles of size c and at a distance r from the initial dipole after evolution through some 
rapidity y. Y is the maximum rapidity determined by the cms energy. The expression for 
the dipole-dipole interaction is given in jSl 13 E] , where it is written as the square of the 
two-dimensional dipole-dipole potential 



,.2 



/(r, c, = -f 



n : , \t + c/2-c'/2\\t-c/2 + c'/2\ 



log 



r + c/2 + c'/2| |r - c/2 - c'/2| 



(36) 



Onium-onium scattering can be formally expressed in the operator formalism presented 
in 0. To that end we define the operators A^(r, b) and A(t, b), where the first one creates 
a dipole of size r and with impact parameter b while the second one destroys the same 
dipole. One usually uses .4.' and A for the right moving onium and X>t and T> for the left 
moving onium. In the operator formalism there are two vertices, the first one, Vi, is the 
vertex for a dipole splitting process while the second one, V 2 is the vertex containing the 
virtual corrections. These are given by 

Vj = a / d 2 r ij d 2 r in d 2 r jn d 2 b ij 5(r ij + r in + r jn ) 2 K \ x 

J in jn 

11 

and 

V 2 = -2a / cPrijtfbykifa )A^r i:j , b^-)^(r y , by) (38) 

where, in (|37j) . we have defined a = a/2Tr. 
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Equation (|37)l is just the vertex for the splitting of dipole ij into in and jn, and here 
the delta function guarantees that the dipoles are connected. The second vertex simply 
takes care of the disappearing mother dipole, which is indicadet by the minus sign. This 
must be taken into account just as we did in section 2.1.1. We define Vr = Vi + V 2 for the 
right moving onium and Vl is defined similarly. The amplitude for multiple scatterings is 
then given by 

F (k \b,Y) = (0\e Al+Vl ^^e yV ^ Y - y ^VUb,r 2 )AUo,r 1 )\0) (39) 

k\ 

where we always choose our coordinate system such that n has impact parameter 0. Here 
|0) is the state defined by .4|0) = £>|0) = 0. Ai is defined by 

Ai= f d 2 bd 2 cA(b,c) (40) 

with a similar expression for T>\. In (I39j) the factor / describes the interaction, the factor 
involving V describes the evolution of the onia states where we create any given number 
of dipoles for all possible positions and sizes. T>^ and A 1 simply create the initial dipoles 
while e" Al+Vl searches for any given number of dipoles for all possible positions and sizes. 
The quantity /, in (|39|h is given by 

f=f d 2 rd 2 r'd 2 cd 2 c'f(r - r', c, 0')!%', c')V(r', c')A\r, c)A(r, c) (41) 

The meaning of the integral is that we count the number of dipoles for all possible positions 
and sizes in the right and left moving onia (T>^T> and A' A are approximately the number 
operators). We then multiply by the interaction amplitude for the dipoles and finally we 
sum all contributions. Thus / has the meaning of a total amplitude. Summing pomeron 
exchanges to all orders one arrives at the formula for the S matrix (which is our goal) 

S(ri,r 2 ,b,Y) = 1 + F(r l7 r 2 ,b,y) 

= (0|e^ 1+I,1 e- / e^ +(y - y)v ^ t (b,r 2 )^ t (0,r 1 )|0) (42) 

To be used in a MC simulation however, this expression for the S matrix is not so good. 
To find a more suitable expression we first note that the probability to find n dipoles in an 
onium state with initial dipole given by &o and r , and at a rapidity Y, is in the operator 
formalism given by 

P n (bo,ro,y) = He yv ^t(b ,ro)|0> (43) 

Assume that we connect the two onium states at the center, i.e y = Y/2. Then the S 
matrix can be written in the following form 

S( ri ,r 2 ,b,Y) = Y,Pa(r,0,Y/2)P a ,(r',b',Y/2)e-^' (44) 
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where the sum runs over all possible configurations, a and a', and f a ^ is given by 

Na N a , 

/«*' = ££/(*<-*$■»*,<) (45) 

i=l j=l 

N a {N a /) is the total number of dipoles for the configuration a(a'). In arriving at the formula 
for the S matrix, it is important that one sums the contributions from the multiple scat- 
terings for a particular configuration before the averaging over all possible configurations 
is done. The reason to this is that the multiple scattering series does not converge, which 
is due to the fact that rare configurations give large contributions to the series, especially 
for higher orders. It is easy to see that the S matrix, as given in (|44j) . satisfies the unitarity 
constraints 

S4;Y, P * P *'< 1 ( 46 ) 

a,a' 

since f a>a > ^ 0. The unitarised amplitude is thus given by 

\F\ = 1 - S = 1 - J2 P a P a' e ~ fa ' a ' = Ys P - P - ,{l ~ e " /a ' a ') ( 4? ) 

a,a' a,a' 

A multiple scattering series for the amplitude emerges when one expands the e~? term 
in the S matrix. Therefore the one pomeron amplitude is given by 

|F«| = 1 -J2 P aPa'(l ~ fa,a>) = £P«P«'(1 " (1 ~ fa,a>)) = ^ P a P a > fa,a> (48) 

a,a' a, a' a,a' 

The factor 1 — e~* a < a ' , in (|47|l . has the meaning of the total reaction probability for the 
onium-onium collision. Therefore all we need to do is to calculate 1 — e~* a > a ' for randomly 
chosen initial conditions, and the unitarised cross-section for onium-onium scattering will 
be given by 

a = / d 2 b(l - e _/ ».«') (49) 



To evaluate the integral one must choose a finite area over which the random positions 
are chosen. This is quite easy to do, since the amplitude decreases very fast for large b. 
By checking the amplitude for large b one easily finds a good value for the maximal value 
of b. In order to obtain the one pomeron amplitude we just do the same thing, but with 
1 — e~f a ' a ' replaced by f a , a ', which is easily seen from (I48j) . Therefore the one pomeron 
cross section is given by 

<X (1) = [d 2 bfa,a> (50) 



Some comments on formula ()36j) are in order here. The four terms in the logarithm are 
the four distances between the ends of the dipoles, see figure 7. If the size of one of the 
dipoles is much smaller than the other n — > r 2 and r 4 — > r 3 and thus the amplitude will 
tend to zero. Therefore small dipoles interact weakly. Also, if the size of one of the dipoles 
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Figure 7: With the notations in the figure the logarithm term in the interaction is log ^ 
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Figure 8: The average number of dipoles with sizes greater than or equal to the initial 
dipole. Although the number of dipoles increase exponentially, the number of large dipoles 
increase very slowly. 
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Figure 9: The average dipole size divided by the initial dipole size as a function of the 
rapidity, y. We see that the dipole size drops rapidly as y increases, just as we expect. 

tends to infinity and one of its ends is not very close to the other dipole, the distances will 
again be mutually equal and the interaction very weak. 

The number of dipoles in an onium state increases exponentially with y. However, as 
y increases it is easy to see that the average dipole size will decrease very fast. This is so 
because p = e~ y /p + gets very small for large y, and the probability to get small dipoles 
increase as p decrease. Since p drops exponentially we also expect the average dipole 
size to drop rapidly. This can be seen in figure El As the number of dipoles increases, 
the interaction term between the two onium states will of course increase, and since the 
number of dipoles increases exponentially one would expect that the cross section also 
increases exponentially with y. However, because of what we just said above, and the fact 
that interactions involving small dipoles are very weak, the growth rate is reduced. In fact, 
the number of large dipoles rises very slowly (with large we here mean dipoles which are 
larger than the initial dipole), see figure |Hl In figure El we plot the unitarised amplitude 
as a function of the impact parameter bo, at a rapidity Y = 14, and for r\ = r2 = r$. We 
see that the amplitude is well below 1 for all impact parameters and that it drops quite 
fast with increasing 6 - We also plot the one pomeron amplitude as a function of b and 
for Y = 14, see figure ^2 It is interesting to see that the one pomeron amplitude is below 
1 for all bo and that it is very similar to the unitarised amplitude. 

This is in sharp contrast to previous studies, for example jH], in which energy conser- 
vation is not included. Here the difference between the unitarised and the one pomeron 
amplitudes is quite large. In |H] one can also see that the one pomeron amplitude is above 
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Figure 10: The unitarised amplitude as a function of b /r , r being the size of the two 
initial dipoles, for Y=l\. 
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Figure 11: The One Pomeron amplitude as a function ofb /r for Y=l\. 
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1 for b < r and that both amplitudes are above zero for 6 ^ 10ro. 

4.2 Onium-Nucleus Scattering 

As a second application we study scattering off a nucleus. As in section 2 we consider a 
nucleus of cylinder shape with some radius B. In the transverse plane the nucleus will look 
like a circle with radius B. The center of the nucleus is located at b . The position of our 
initial dipole is kept fixed, its coordinates are given by r = (— 0.5ri , 0) and ri = (0.5ri , 0). 
The impact parameter is then simply equal to bo- In order to calculate the cross section 
for this event we have to find an expression for the interaction between the dipoles and the 
nucleus. First of all one should note that the amplitude, for dipole-nucleus scattering, is 
going to depend on five parameters. Two of the parameters give the position of the center 
of the nucleus in the transverse space, namely &o = |b | and ijjf,, the third parameter is 
the size of the initial dipole, ro, the fourth is the angle describing the relative orientation 
between the target and the initial dipole and finally the fifth parameter is the rapidity. 
The relative orientation however, is not needed when the nucleus is cylindrical, which is 
quite obvious. The rapidity and the size of the initial dipole come in naturally in a MC 
simulation and are not needed explicitly in the interaction formulas. Therefore we are left 
with two parameters instead of five. 

As mentioned above we picture the nucleus as a circle in transverse space. We will 
view the nucleus as a collection of color dipoles and the interaction between an incoming 
dipole and the nucleus will be given by the interaction between the incoming dipole and the 
target dipole in the nucleus, convoluted with the dipole distribution of the nucleus. The 
interaction formula between two dipoles is again given by (J36|) . The dipole distribution can 
be chosen in a variety of ways but not arbitrarily. First of all, we want the dipole density 
to go to zero as the dipole size goes to zero. This is because, as mentioned in section 3.1, 
small dipoles correspond to hard gluons and if we allow too small dipoles in the nucleus 
we will violate energy conservation. The dipole density should also go to zero when the 
dipole size approaches infinity, we do not want dipoles which are larger than the nucleus. 
The dipole density should also be large at the center of the nucleus, an incoming dipole 
with zero impact parameter should have a larger chance to interact with the nucleus than 
a dipole with non zero impact parameter. Let b« be the impact parameter of the target 
dipole i , i.e a dipole belonging the nucleus. The position of dipole i relative to the center of 
the nucleus is then simply given by b = b^ — b . We also expect the dipole density to drop 
quite fast for large b. With these considerations we have chosen the following distribution 

dV = d\d 2 br^e- r 'e-^ (51) 

Note that this whole expression is dimensionless, the r, and b are all divided by a unit 
parameter with the dimension length. Of course, in this expression, there should also be 
a normalization constant, which may come from energy conservation. We are here less 
interested in the precise value of the amplitude however, primarily we want to investigate 
how the amplitude behaves. In future studies a normalization constant should be included. 

24 



In our calculations we have studied an example where B is equal to five unit lengths. We 
see that the nucleus will consist of a collection of dipoles, with sizes small compared to the 
nucleus, and a higher density near the center. The interaction between an incoming dipole 
j and the nucleus is then given by 



b - b , r i5 Tj) = / dVf(hj - b - b , r u r,) 



/( bj - b - b , Ti, r d ) = I dVfihj - b - b , Ti, r d ) (52) 

Once we have this expression we just sum over all incoming dipoles j and the total 
cross section is, in analogy to ([49)1 . given by 



a 



/ d 2 b l-expl-^2 f(bj - b - b , r i} r,) J 



(53) 



The amplitude contains multiple pomeron exchanges up to all orders and it satisfies the 
unitarity bounds. We calculate (|52j) and (J53)l using MC techniques. For this one has to 
choose a finite area in which the last integral is evaluated. This is quite easy to do since the 
amplitude falls off quite fast for large impact parameters. It is usually more than enough 

to take b 0max «5(r + 5). 

4.3 Results 

In this section we present the results we obtained from our program. To start we like 
to point out that all calculations have been performed using a fixed coupling constant, 
a = 0.2. Although a running coupling constant would obviously be more justified, we have 
avoided it in order to make the whole Monte Carlo program simpler. As we described in 
the last section, the main objective has been to obtain the cross-sections for the events we 
studied. Therefore we start by presenting the cross-section for onium-onium scattering for 
different configurations. 

In figure El we compare cross-sections obtained from one-pomeron and unitarised am- 
plitudes. We see that the difference, hence the saturation, is quite small. We could observe 
this already in section 4.1, where we compared the amplitudes as functions of the impact 
parameter. Observe that the scale in figure El is not logarithmic like in other figures. 
If plotted on logarithmic scale the difference between the cross-sections would barely be 
visible. One can see that the difference increases with Y, this is no surprise since / will 
increase with Y and therefore the difference between / and (l-exp(-f)) will also increase. 

Figure ED shows the cross-sections for different initial conditions obtained from our MC 
program. The difference compared to previous studies, [8] , is clear, our cross-sections grow 
much more slowly. One of the reasons for this is, as explained in section 4.1, that the 
average dipole size decrease rapidly as Y increases and since small dipoles interact weakly, 
the growth of the cross section is reduced. The other reason is discussed below. 

In figure El we show a/r 2 as a function of the rapidity calculated for onium evolution 
with constant ultraviolet cut-off p, i.e without energy conservation. Comparing with figure 
El we clearly see that a grows much faster in this case. 
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Figure 12: Comparing the cross-section calculated from one pomeron and unitarised am- 
plitudes, as a function of the rapidity Y . 
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Figure 13: ojr\ for various initial conditions as function of the rapidity Y. Here r\ is kept 
fixed while we vary ri- All calculations are done using the unitarised amplitude and for 

y = Y/2. 



26 






1 1 
r Constant p=0.02r + 

+ 

+ 
+ 


i 

+ 


+ 


i 

+ 


+ 


- 


1 

1 


+ 
+ 


i 




i 




: 



10 



Y 



12 



14 



Figure 14: o~jr\, r$ being the size of the two initial dipoles, using the unitarised amplitude, 
as a function of rapidity Y and without energy conservation. Note that the cross-section 
grows more rapidly compared to figure 13. 
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Figure 15: A, plotted as a function ofr 2 /ri, calculated from the unitarised cross-section for 
onium-onium collisions. 

It is again interesting to see that the unitarised cross-section in figure El grows even 
faster than the one pomeron cross-sections in figure El we could see a hint of this already 
in section 4 where we plotted the amplitudes. We see that energy conservation has a non- 
negligible effect on the growth of the cross section. It seem as if one does not even has to 
use the unitarised amplitude, at least for the Y intervals we investigated. 

Recall from section 2.1.2 that the gluon distribution function in the BFKL region is 



given by G 



x 



v . Here lnl/x ~ y, which gives G ~ e y . Therefore the BFKL equation 



predicts an exponential rise of the cross section, a ~ e Xy . Thus lncr ~ Xy where, in the 
leading log approximation, A is given by a " A ln2= 4aln2^ 0.55 for a = 0.2. By trying 
to fit a straight line for lna between Y — 10 and Y = 19 we have obtained values for A. 
The results are shown in figure ED We see that the A we have obtained is much lower 
than Xbfkl- We explained one of the reasons for this above. The second reason is that 
the steps in y, Ay, gets larger when we have energy conservation, which results in fewer 
dipoles. That Ay gets larger can be seen from the phase space diagram in figure El Since 
we have both p + and p~ conservation the available phase space for emitting a gluon with a 
small Ay is small. Therefore it is more likely to take bigger steps in y. If Ay is the typical 
rapidity interval where the number of dipoles double, we can estimate the total number of 

dipoles with iV ~ 2 Y ^ Ay . This can be written as iV 



exp 



^-Y ) . Therefore we get A 



hi2/Ay. It is not really true that a increase proportional to iV as the average dipole size 
is also important. The large Ay does, however, give a very essential contribution to the 
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Figure 16: o jr\ for unitarised and one pomeron calculations for onium-nucleus scattering 
as function of the rapidity Y, for B/ro = 10. 

reduction of the parameter A. Thus it is expected that we obtain a lower A than in jHJ. 
This is also expected from theoretical considerations. Taking into account higher order 
corrections reduces the value of Xbfkl- For values, as ~ 0.2, the value for A in NLO is 
even negative. This is of course not so realistic and there are different calculations that 
give some other values but it gives an idea of how A behaves. A large fraction of these 
higher order corrections are related to energy conservation, and therefore it is expected 
that one obtains lower values for A. Observe that our values for A are not exact, there 
is some uncertainty in trying to fit a straight line for lna. We have calculated A also for 
constant p, and, for p = 0.02r , r being the size of the two initial dipoles, we get A = 0.32. 
This is a bit lower than Xbfkl = 0.55 but the value of Xbfkl is valid in the limit p — > 0. 
Figures ITBI and IT7I show the results for onium-nucleus scattering. The saturation effects 
are more visible here. For small Y we see the power like rise of a and as Y increases the 
cross section starts to saturate. The reason that we can see the difference between the 
one pomeron and the unitarised amplitudes even for small Y is because the onium now 
interacts with a much denser object. Observe also that the Y in onium-onium collisions 
and the Y in onium-nucleus scattering has different meaning for the evolution of an onium 
state. In onium-onium scattering events, the onium states are evolved up to Y/2 while 
in onium-nucleus scattering, the state is evolved up to Y. For a low rapidity, such as 
Y = 2 (Y = 4 in onium-onium events) the onium state often consists only of the initial 
dipole, therefore an onium-onium collision will just mean a collision between two dipoles. 
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Figure 17: o jr\ for unitarised and one pomeron calculations for onium-nucleus scattering 
as function of the rapidity Y, for B/r® = 5. 

Hence it will be very difficult to see any difference between the one pomeron and unitarised 
amplitudes. In nucleus collisions however, that single dipole will collide with a much denser 
object, which consists of many dipoles, and therefore the amplitudes will not be so small 
and the difference will be more visible. There is a difference here compared to the case of 
constant p. If one uses a constant p , which is sufficiently small of course, the steps in y 
will be smaller and therefore there will always be additional dipoles even for smaller Y. 
Hence the difference in the amplitudes will be visible, though small. This can be seen in 

5 More on Saturation 

5.1 The Dipole Fusion Factor 

In this section we will focus on a process that we excluded in our analysis, namely the 
question of dipole fusion processes, the creation of a dipole by destroying two. To make 
the analysis easier we abandon energy conservation and switch to a constant cut-off, p. 
We mentioned the fusion process a little when we presented the Balitsky-Kovchegov equa- 
tion, (|21|). As we said, the BK equation does not contain this possibility. This is so because 
Mueller's Dipole formulation does not contain it, and the BK equation is obtained from (J2U|) 
by summing multiple pomeron exchanges up to all orders. The inclusion of dipole fusion 
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Figure 18: A dipole fusion process. The dipoles denoted by c and d go together and form 
T\ and the evolution continues 



processes in the wave function formalism is not an easy task and is yet to be solved. 

Let us return to Mueller's formulation. As we recall from section 2, Mueller's model 
gives us the squared wave function for an onium state with arbitrarily many gluons. We 
also saw that, when two new dipoles, in and jn, were created as the result of a gluon 
emission from dipole ij, there was a factor adrx n 2 % \ associated with this splitting. The 

X in X jn 

wave function, 4>^ n \ is then obtained from (f)^ ' by considering all possible ways in which 
we can arrive at a state with n gluons. 

Look at figure QBl where we follow some evolution that leads to the formation of the 
dipole labeled by r . Then r splits into a and b , where a in turn splits into / and c, and b 
splits into d and e. Then, from what we said above, we will get a factor (denote the gluon 
at the center with A) 

(f rA I^-— — (54) 

a?b 2 f 2 c 2 d 2 e 2 

for these splittings (of course we will also get factors, d 2 r, from the gluons on the top right 
and bottom right corners, but these factors are not interesting for our analysis). Note that 
in the large N c limit only dipoles that share a common gluon, i.e neighboring dipoles, can 
fuse. After these splittings we let c and d fuse and form r\. The fusion process will be 
associated with some weight Z ■ dy, and the question is what this factor looks like. Observe 
also that we do not need to write r 2 , in (|54j) since there will be a factor 1/Vq from the time 
ro was formed and therefore these factors cancel. Instead we can include a factor r \ which 
will be present when the dipole t\ emits a gluon, later in the process. Thus we get 



d 2 r A 



f 2 c 2 d 2 e 2 



Z(n;c,d)r 2 



(55) 



where the arguments of Z indicate that the dipole r\ is formed out of the dipoles c and d. 
As discussed above, we require that the theory looks the same if we view the process from 
the other end. If we start from the right and follow the evolution we will see r\ form first, 
then 7*1 will split into d and c, c will in turn split into a and / while d splits into b and e, 
and finally a and b go together to create r and the process continues. Now we will get the 
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factor 

d2rA ^W¥^ z{r ^ h) (56) 

Applying the same arguments as above we get rid of the factor r\ and we can include 
r\. Hence we obtain 

rf2r ^7^5v z ( r o; a ' 6 ) r o ( 57 ) 

By the symmetry argument (|55|) and (foTj) should be equal, hence 

r o err.. . „ ^ r ? 



Z(r ;a,b) = Z(n;c,d) 



a 2 b 2 e 2 f 2 ' ' c 2 d 2 e 2 f 

^Z(r ;a,b) = -^-Z(r x ;c,d) (51 



o?b 2 ' c 2 d 2 

Most cascade models are semi-classical but despite this they are very successful in 
describing experimental results, such as in e + e~ annihilation. The DGLAP region in DIS 
is also successfully described despite the fact that DGLAP, which contains probabilities 
instead of amplitudes, is not fully quantum mechanical. In this, semi-classical, spirit we 
assume that the fusion factor, just like the splitting factor, is local and only depends 
on the dipoles involved in the fusion. A full quantum mechanical treatment would give 
contributions which include interference factors but, given the success of these semi-classical 
models, we will primarily try a locally factorizing approximation. The assumption that 
the fusion factor is local means that Z(r§\ a, b) (Z(r\\ c, d)) only contains the lengths r , a 
and b (r 1 ,c and d). Thus the only possibility gives 

a 2 b 2 c 2 d 2 

Z(r ; a, b) oc — — and Z{r\ ; c, d) oc — ^- (59) 

2 2 

Generally, if dipoles r in and r jn go together and form r^ we get a factor "Y" . To make 

the formulas in (|59j) complete we note that the fusion factor Z must be dimensionless. 
Therefore the expressions in (|59|) should be divided by a factor which has the dimension 
of an area. A reasonable choice would be to take the area of the triangle that is formed 
by the two disappearing dipoles and the dipole that is created. If this is done however, 
we get a factor 1/A\ on the LHS of (|55|) while we get I/A2 on the RHS, where A± is the 
area of triangle formed by a, b and ro while Ai is the area of the triangle r\cd. Since in 
general A\ 7^ Ai, the equality, hence the symmetry, will be violated. In order to keep 
the symmetry one can easily see that we must have something that looks like \j\J A\Ai- 
Including A\ {A%) in the fusion factor, Z[v\\ c, d) (Z(tq\ a, b)), means that we remember the 
steps that we took before. From what we said above however, we avoid this possibility and 
assume that the fusion factor is local. Therefore Z{j\\ c, d) (Z(tq\ a, b)) will not contain A\ 
(.4.2)5 but as we just saw we cannot get a theory which is symmetric by only including A2 
(Ai). Thus the only possibility is that the fusion factor contains a parameter, which has 
the dimension (length) 2 , and is a fundamental parameter in the theory. One can ask what 
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Figure 19: On the left side we see b and c fuse to create d at a rapidity y and then a splits 
into e and f at y' . The figure on the right side shows the same process viewed backward. 
Now e and f undergo fusion to create a at y' and d splits into b and c at a rapidity y. 

kind of fundamental scales there are. The fundamental scale of QCD is Aqcd which has 
the dimension of energy. Therefore one can make the heuristic argument that Z oc A.q CD . 
It is also reasonable to assume that the fusion process, just like the splitting process, is 
associated with a factor a. The complete formula for the fusion of dipoles in and jn into 
ij is then given by 



^yiji'ini ' jnj 



T; 



«£ 



in' jn A 2 

2 ^QCD 



T 



(60) 



ij 



where £ is a dimensionless, free parameter in the theory. 

5.2 Putting the Fusion Factor in a MC Program 

Now that we have the fusion factor we can use it to generate new y values, which we denote 
by y. As mentioned in section 3 we generate, for each existing dipole, a y value and then 
we accept the gluon with the lowest ?/-value. For each dipole we also generate y-values 
for possible dipole fusions and if the lowest y happens to be lower than all the other y we 
accept the corresponding fusion process. The y-values are given by 



y = y 



a £ Aqcd ■ *»■ j,> 



— ^ — ^-lni?, 



(61) 



where R is a random number and y is the maximum rapidity before the fusion process. 

One may ask what physical meaning y has. Consider the processes on figure El where 
we picture the evolution of an onium state for some rapidity intervals and then we view 
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the same steps backward. When we come to y we generate and choose y and let b and c 
go together to form d. When viewed backward the meaning of y will be the rapidity at 
which d splits into b and c. The same thing goes for the process where a splits into e and 
/ at y' . Viewing the process backward we see e and / undergo fusion and create a and 
the rapidity y' would then be the rapidity which we would have obtained from (jHTj) for the 
backward evolution. 

5.3 Stability 

Let us return to figure HH1 and the gluon at the center denoted by A. What happens if 
we vary r^? In particular consider events where A is emitted far away from the mother 
dipole. After the fusion of c and d into r\ there will be no terms left that depend on r^. 
Therefore, if we integrate d 2 rA over large distances we will get an infinite contribution. The 
problem is that even though the probability to emit the gluon A at very large distances 
goes to zero, the probability of having the fusion event c + d — > r\ goes to infinity. This 
can be easily seen from the fusion factor S-f-. Hence in order to make the theory stable the 
contribution from very large dipoles must be suppressed. Enforcing energy conservation we 
got a cut-off for the large dipoles from p_ conservation but in the case of constant p there 
is no limit for how large a dipole can be. It should be reasonable to use an infrared cut-off 
even if one does not consider fusion processes. This has to do with the Froissart bound, 
which says that the cross section cannot grow faster than ln 2 s. The Froissart bound follows 
from the assumption that a field, which acts with the exchange of massive particles, has a 
finite range. Therefore one should cut off contributions that come from too large distances, 
which means that it is not right to allow too large dipoles. In QCD this is expressed in 
the confinement mechanism, which suppresses large dipoles. A reasonable infrared cutoff 
is given by 1/Aqcd- in a Monte Carlo program the easiest way to implement this would be 
to generate the y, r x and r y as before and afterwards accept these values with a probability, 
V ~ e -( r i+ r 2)A<2c.D_ Here r\ and r2 are the sizes of the dipoles which were created in the 
splitting process. 

Next, let us study what happens when the gluon A is emitted very close to one of its 
mother partons. Denote the smaller distance to the mother partons by 5. Then for very 
small 5, the probability to emit A gets very large but in order to have a fusion process, 
involving a dipole that contains A, we must first have an emission from the dipole with 
length S, which then will get a weight proportional to 5 2 . Therefore the total probability 
for the whole process will not be large, and since the value of 8 varies over a small space 
we will not get a divergent contribution. Therefore we see that including fusion processes 
causes no unpleasant instabilities in the theory. 

One can ask how large effect the fusion processes will have on the onium evolution. 
Since the probability to create small dipoles is high (of course we assume that p is small) 
there will be a lot of small dipoles. On the other hand, the probability to have fusion 
processes involving small dipoles is very low and therefore it would seem as if the fusion 
effects would not be so visible. Small dipoles interact very weakly however, therefore it 
does not make so much difference if many of them disappear due to fusion processes or 
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Figure 20: o jr\ plotted as a function of Y for various £ and for constant p = 0.02r . 

not. The largest contribution to interactions come from the larger dipoles and these have 
non- vanishing probabilities to undergo fusion, therefore for a state with many large dipoles 
the effects of fusion processes should be visible. 

5.4 Results 

In this section we present the results we obtained from our MC simulations including dipole 
fusion processes. We have run the program for different £ and studied how the cross section 
behaves for these £- values, in onium-onium collisions. In figure I27H we see some of the result 
obtained by including fusion effects. We choose r such that tqKqcd = 0.166. The cross 
sections were calculated for fix p. We see that A decreases as £ increases and lncr still seems 
to depend linearly on Y . It doesn't seem as if the growth starts to saturate but of course it 
might be that the saturation sets in for larger Y. This is actually expected in the large N c 
limit where only neighboring dipoles can fuse. From these result we have also calculated 
the different A and in figure l2"Tl A is plotted as a function of £. 

Let Xf denote the A-value we obtain from onium evolution including fusion processes. 
Similarly we introduce \ec, where EC stands for energy conservation. From figure I2~T1 we 
see that A^ drops slowly for small £. This is not strange since for small £ (£ < 1), any 
fusion process will have very a small probability to occur, and changing £ will not affect this 
much. As £ increases we see that A^ starts drop more sharply since the fusion processes 
become more probable. We have not gone any further than £ = 1000 since running the 
program for larger £ is extremely time consuming. The reason is that the steps in y gets 
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Figure 21: A plotted as a function of £. 

extremely small when £ is very large, which follows immediately from (jfilj) . 

Looking at the figures E3 and ED we see that A^ > Xec- This means that energy 
conservation has more effect on the growth of a than fusion processes. It can be that the 
fusion processes become more important than energy conservation at higher rapidities. Of 
course when estimating Xec we used values within the interval 10 ^ Y ^ 19 but for Xp 
we used values between Y = 8 and Y = 14. The results would be slightly modified if one 
uses this interval to determine Xec- F° r the lower Y range we would, for r± = r 2 , find 
Xec — 0.12 instead of Xec — 0.10. The new Xec still satisfies Xec < A^ though. 

If we calculate the cross section without fusion processes but imposing the infrared cut 
off we get A ~ 0.21. This might seem strange since A^ for £ = 0.1 is approximately equal 
to 0.22. This is not surprising though since there are some uncertainties in the A-values. 
Also, Xf((, = 0.1) should more or less be equal to A, since any fusion process will be highly 
unlikely for these values. Especially when we have an infrared cut off. Finally we note that 
the fusion processes have little impact on A for £ < 10. 

6 Outlook 



In this section we will present three problems that can be investigated in future studies. 
The first problem is certainly something very interesting. Consider two onia states that 
we connect at some rapidity, see figure 1221 Denote the right (left) moving onium with & 
(=£?). Assume that the interaction between the two states occurs between the dipoles a 
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Figure 22: The two states are connected between dipoles a and b. The dots just indicate 
that there has been an evolution up to a and b and & (££) denote the right (left) moving 



onium. 






Figure 23: In this figure, c and d are the diagonal distances between a and b. 

and b, where a (b) belong to 3% («£?)• The formula we used for the interaction between two 
dipoles is given by l)flfij) . This formula can be written in the form 



a,. 2 cd 
— hr 



<r, c 2 d 2 



2 

'In 2 ' 



ef 8 e 2 f 



(62) 



For the meaning of c, d, e and /, see figure|2Sl Observe that figures 22 and 23 are extremely 
simplified, just to give an idea of what happens. In general it is highly unlikely that two 
onia collide head on and the dipoles a and b can have any relative orientation. 

The expression above does not factorize in the same way as the rest of the splitting 
factors in the chain. However, if the separation between the dipoles is large compared to 
their sizes it does reproduce an expression that is similar to the splittings factors. Thus 
e and / are much larger than a and b. We can choose the vectors such that c = a + e, 
d = e + b and f = a + e + b. Then we have 



c 

d 2 

2 



./ 



a 2 + e 2 + 2ae cos 9 
b 2 + e 2 + 2be cos <\> 
a 2 + b 2 + e 2 + 2ae cos 9 + 2be cos + 2ab cos ip 



(63) 



where, using our definitions of the vectors, we have ip = 9 — <f>. Of course one can choose 
the orientations such that ip = 9 + <fi also, but this will not change the result. Next, we 
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define e — - and S — -. Then we get 

p. p. o 



z 2 
e" 2 

e 



: 1 + e 2 + 2e cos 6 1 



2 ^ 1 + 5 2 + 25 cos 



c2 

,2 , x2 



^r - l + e 2 + 5 2 + 2ecosfl + 2£cos0 + 2e<5cos^ (64) 

e 2 



Using these relations we obtain 



,2 



2 , o,„„„/Q ± /,2 , o,„„„/Q\2 



In— w e +2ecos0 (e 2 + 2ecos( 

r/ 2 1 

In— « 5 2 + 2<5cos0--(5 2 + 25cos0) 2 
e 2 2 



f 2 
ln^- « e 2 + 5 2 + 2ecos6> + 25cos0 + 2e<5cos^ 



-.-. r 2 ' x2 
1 



--(e 2 + 5 2 + 2e cos + 25 cos + 2e5 cos ^) 2 (65) 



Therefore we get 



c 2 d 2 c 2 d 2 f 2 

ln^r— = In— + In— - In— 
e j e e e 

Ki -2e5 cos ip + 2t 2 5 2 cos ty + e 2 5 2 + 2e 2 5 cos + 2e 3 <5 cos tp + 2e5 3 cos ip + 

+ 2e5 2 cos + Ae5 cos cos # + Ae 2 5 cos # cos -0 + 4e<5 2 cos ?/> cos 

~ 4e5cos0cos# — 2e5cos-0 (66) 

where we have kept terms to second order in e and 5. Using the relation ip — $ — we get 

c 2 d 2 
In sa 2e£(cos0cosfl — sin sin 0) (67) 

e 2 / 2 

Squaring this expression we get 

Ae 2 5 2 (cos 2 cos 2 9 — 2 cos sin 6? cos sin + sin 2 sin 2 0) (68) 

If we consider all the possible orientations of a and b, while keeping e fixed, we should 
take the average of the expression above 

4e 2 5 2 (cos 2 cos 2 9 - 2 cos sin 9 cos sin 9 + sin 2 sin 2 9) 9>(j) = Ae 2 5 2 ( - + - ) (69) 



Hence we have 

8 ~"e 2 / 2 8"~ " V 2 ) < 2 P " r-'f 
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a 2 c 2 d 2 a 2 fa \ 2 a 2 b 2 a 2 b 2 



Thus with the approximation that the separation bewteen the interacting dipoles is much 
larger than their sizes we can write the dipole-dipole amplitude as 

a 2 b 2 
-c?% ( 71 ) 

Let us return to the onia states in figures E21 and 1221 It is easy to see that, using (|52j) . we 
do not get the same result if we connect the states at different rapidities. This is somewhat 
disturbing, we would like to have an amplitude such that we can connect the onia states at 
any given rapidity and get the same results. This would be true if the total weight for the 
evolution of the states £% and ££ , multiplied with the amplitude between a and b would be 
equal to the weight of a single chain, which is spanned between the ends. The weight for a 
chain is simply the product of all the splitting factors (to make things simpler we assume 
no fusion processes). To make the weight symmetric the reader can easily verify that we 
must have an amplitude that is equal to (|7Tj). 

It would be interesting to find out what the splitting factors would be, if one wants to 
have the same symmetry by using the amplitude (|62[). Obviously such a splitting factor 
would be similar to (|62|). It would perhaps follow from next to leading order calculations, 
which is a better approximation of Mueller's model. Then, just as we did in deriving 1)71)1. 
it would be possible that one obtains the splitting factor in Mueller's model as an approxi- 
mation of the new splitting factor. In this new splitting factor, the contributions from the 
different dipoles would not factorize as before, but it would contain different interference 
terms, which can be seen from (162)1 . 

We leave the first problem and go on to the second. As mentioned earlier we have tried 
to construct a theory that is right-left symmetric. The weights we use are symmetric but 
everything is not exactly symmetric. Consider figure |^ where the rapidity increases to 
the right. Recall from section 3.2 that the transverse momentum of a gluon is decided by 
the shortest distance to another gluon, with which it has formed a dipole. 

We look closer at the gluon denoted by a, as we can see this gluon has different links 
connecting it to the other gluons in the two figures. Therefore it can have different p± 
values depending on which side we start from. When we evolve from left to right the 
dipole denoted by b forms at some stage but when we start from right, the same dipole 
never forms, instead we get the dipole c which is not seen on the left figure. Observe that 
the weights will be the same in both the right and the left figures since there will be no 
factors depending on b and c. Therefore it is not a big problem that we do not have the 
Pj_ symmetry. An improvement would be to make a better approximation, in choosing the 
Pj_ values for an emission, that respects p_i_ conservation. Another, and a very natural, 
improvement would of course be to use a running coupling constant a, which can be an 
interesting project for the future. 

Finally, we look at the third problem which might also be interesting. Going back to 
the BK equation, (|2~T)l . we mentioned that this equation was derived from ()2*U)l which uses 
a constant ultraviolet cut-off p. It would be interesting to try to modify the BK equation 
with a energy conserving, i.e y dependent, cut-off. One obvious change would be in the 
Sudakov factors of (|20[). As we mentioned in section 2.2 the first Sudakov factor of ()20|) 
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Figure 24: On the left figure we start from the left and evolve towards the right while on 
the right figure we the evolution starts at the right end and stops on the left 



described the probability that nothing happens during the evolution, from y 
This Sudakov factor has the form 
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to y = Y. 



(72) 



for constant p. If we make the cut-off y dependent we get 
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(73) 



where, recall from section 3.1, p + is the positive light cone-momentum of the onium. So the 
Y dependence of the Sudakov factors changes slightly. The second factor in (J2*U|) contained 
the probability that nothing happened after the last gluon is emitted at a rapidity y, which 
can be anywhere between and Y. This Sudakov factor is then given by, for a constant p 



S = exp 



Y 



—a I dy 1 d T2 



j2„ r io 



'12' 20 



exp 



-2oln^(y-y) 



(74) 



Now one could think that we just do the same as we did above when we switch to a y 
dependent cut-off. That is true but there is also a problem, we must know what p + is, 
now it won't be just given by the energy of the onium. This problem will also appear in 
the X2 integral in (|2*U|) . p(y) will now change in every emission and for every dipole we will 
have a different p. Therefore one would need to keep track of all the emissions during the 
evolution, this is not a problem in a MC simulation, after all that is what we have been 
doing in this thesis, but for a theoretical approach it is a major problem (for a constant 
p, these problems does not appear since it doesn't matter how many emissions there have 
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been or which dipole is about to emit, the cut-off is the same anyway). It seems to be 
difficult to write down an equation for a generating functional as in ()20j) . However, it would 
be interesting if this was done. If one could just find such an equation, deriving a modified 
BK equation would not be difficult. Of course the new equation should not be too difficult 
to analyze, if one finds an equation that is not even numerically solvable it would not be 
so interesting. 

7 Conclusions and Summary 

Let us summarize the thesis. In other studies energy conservation has frequently been 
found to have a quantitatively large effect, and also to correspond to a significant part 
of NLO corrections. Our main objective has been to study saturation effects in DIS. To 
do this, we have constructed a Monte Carlo program which is based on Mueller's dipole 
formulation and respects energy conservation. Using our model we have studied A, which 
is, in the BFKL approximation, defined by the relationship a ~ e Xy , and seen that there 
is a clear difference compared to previous studies, such as |H|, where energy conservation 
is not included. Throughout the thesis we have been working in the large N c limit. 

In onium-onium collisions we have seen that the difference between the one pomeron 
and the unitarised cross sections are quite small. It is hard to say if the cross section 
saturates but its growth with energy is reduced and we obtain A in the interval 0.10-0.15. 
In onium-nucleus collisions the difference between the unitarised and the one pomeron 
cross sections are larger, and the growth of the cross section slows down quite rapidly. It 
is also clearer in onium-nucleus scattering that lna does not depend linearly on Y, except 
for small Y. The saturation sets in quite early and the growth of the cross section becomes 
smaller and smaller. In the onium-nucleus case however, one should be quite careful before 
drawing any conclusions. This process was not studied in detail in this thesis but it is going 
to be studied more deeply in future investigations. As we mentioned in the text the dipole 
distribution for the nucleus should be modified by taking into account normalization. The 
distribution we used was a very simple one and can be improved. 

In the second part of the thesis we studied the effects of dipole fusion processes. Using 
symmetry arguments, and in a factorizing approximation, we proposed an expression for 
the fusion factor, presented in (JoTij) . We used this expression in our Monte Carlo program, 
calculating cross sections for onium-onium collisions. Due to lack of time we only studied 
the effects on the unitarised amplitude and without energy conservation. The value for 
the fusion parameter, £, cannot be determined by theoretical arguments however, and we 
have here considered various values for this parameter. Varying £ we obtained different 
values for A, and we saw that fusion processes are not as important as energy conservation, 
in slowing down the growth of a, at least for Y ^ 14. Before any definite conclusions 
are made though, the program should be tested more and it should be run with better 
statistics. It would also be interesting to study larger rapidity intervals, especially when 
including fusion processes. Fusion processes become interesting when the gluon density is 
high, which occurs at high rapidites. Studying larger Y intervals is very time consuming 
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however, especially for onium evolution without energy conservation. The problem is that 
the number of dipoles becomes very large as Y increases, which limits the F-range available 
for numerical studies. 

We should also repeat that one needs to use a large £- values, larger than 10 at least, 
to see any significant effects of fusion processes. Since we derived (ffiT)|) in the large N c 
limit, which means that only neighboring dipoles can undergo fusion, the effects of fusion 
events are not so large for small and moderate Y, (Y < 14). Energy conservation has a 
larger effect on the growth of the cross section and its effects are seen earlier than fusion 
processes. Therefore energy conservation should have higher priority than fusion processes 
in DIS investigations. 
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A Appendix 

In this appendix we derive the formula (J34J1 . To find the y distribution we use the for- 
mula ()31j) . Remember that r^ = (0, 0) and Tj = (1, 0). Therefore we can write the splitting 
factor as 

r 2 . r 1 

' 91 — / -L 



a / d 2 r n ^\- = a / dr x dr y = / dr x dr y f(r x ,r y ) = F{y) (75) 

J rf n rj n J {.rl + rl){(r x -l) 2 + r 2 y ) J 

The integration is done over all of M 2 except for two circles with radius p(y) where the first 
one is centered at and the second one at (1,0). The rapidity dependence in F comes 
from the fact that p = p(y). Using the veto algorithm the y distribution is given by 

y = F-\F(y t )-\nR) (76) 

Here F is the primitive function of F, yi is the initial value of y and R is a random number. 
Finding F~ l is easy as long as p is sufficiently small. But the integral in ()75|) is difficult 
to evaluate when p is not small and since we have a running p we cannot always be sure 
that p is sufficiently small. Therefore one needs to use an auxiliary function g, satisfying 
g ^ /. Before we do this, we note that f(r x ,r y ) is symmetric around r x = 1/2. Therefore 
we do not have to consider the whole M 2 plane but it is enough to look at the plane where 
r x < 1/2. We choose 

9{TX ' Ty) = (^ + r 2 )(r 2 + r 2 + 0.25) (7?) 
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We easily see that g{r x ,r y ) ^ f{r x ,r y ) for all (r x ,r y ) with r x < 1/2. Then we have 

dr x dr y g(r x ,r y ) = 167raln— = 87ra;ln(l H — — ) (7£ 



We use g(r x , r y ) to generate the r x and r y values, which can be easily verified. Next we 
introduce a second auxiliary function, h(y) given by 



e 2 V e 2y r> 2 

h(y) = 8naln(e 2 y + —j!±) > 87raln(l + —j!±) = g(y) 



(79) 



where the inequality follows from y ^ 0. Finding the primitive, H, of h is easy, and it is 
also easy to find the inverse H~ 1 . Using (|7f)jl. with H, we immediately obtain (|32)1 . To 
obtain the correct distribution we must accept the generated r^,^ and y values with the 
probability 

f(r x ,r y )g(y) 

g(r X} r y ) h(y) 

We are not done however since we considered only the plane where r x < 1/2. When we 
generate the r x , we reflect them around the symmetry axis, r x = 1/2, with a 0.5 probability. 
Therefore instead of g(r x ,r y ) we should in (JKUj) use 



g{r x ,r y ) 
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+ 
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(r* + r 2 y )(r 2 x + r 2 y + 0.25) ((r x - l) 2 + r 2 y )((r x - l) 2 + r 2 + 0.25) 
Using this g, the correct formula is given by 
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which gives 
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